#!/bin/bash
set -e

while getopts  ":i:p:" opts
do
        case  $opts  in
        i) 
                        intervals=$OPTARG
                        
				;;
		p)
                        out_prefix=$OPTARG
                        
                ;;
		\?)

				echo `basename $0` [-i interval -p out_prefix]
				exit 0
				;;
        esac
done
shift $(($OPTIND - 1))

#-----------------------------------------------
#-----------------------------------------------
. /mnt/ilustre/app/medical/tools/.var #---------
#-----------------------------------------------


echo interval file is: $intervals
echo out prefix is: $out_prefix

if [ -z $intervals ]; then
        intervals=${data_path}/intervals/1/brca.qiagen.NGHS-001X-Covered.b37.bed
        echo interval file is: $intervals
fi





###--------------- argument may be changed ---------------###


# genome_name=b37.fa
echo using $genome_name as the reference genome 2>>$log 1>&2
# genome_assembly=b37
echo genome assembly is: $genome_assembly 2>>$log 1>&2
# dbsnp_version=138
echo dbsnp version is: $dbsnp_version 2>>$log 1>&2

# data_thread_num=8
# cpu_thread_num=4

# java_memory=16g
echo java memory: $java_memory 2>>$log 1>&2

# snpeff_db_version=GRCh37.75
echo snpeff database: $snpeff_db_version 2>>$log 1>&2


echo ref genome is: $ref_genome 2>>$log 1>&2



reads_seq_1=$1
reads_seq_2=$2

echo reads are: $reads_seq_1 and $reads_seq_2 2>>$log 1>&2




echo 2>>$log 1>&2
echo 2>>$log 1>&2
echo gatk UnifiedGenotyper 2>>$log 1>&2
java -Xmx$java_memory -jar $gatk \
	-T UnifiedGenotyper \
	-R $ref_genome \
	-I 1.realn.recal.bam \
	-o 1.snp.vcf \
	-glm SNP \
	-baq CALCULATE_AS_NECESSARY \
	--baqGapOpenPenalty 40.0 \
	--downsample_to_coverage 2000 \
	-nt 8 \
	-nct 2 \
	-stand_call_conf 50.0 \
	-stand_emit_conf 30.0 \
	--min_base_quality_score 20 \
	-contamination 0.0005 \
	--max_deletion_fraction 2 \
	-L ${intervals}
	
	
echo 2>>$log 1>&2
echo 2>>$log 1>&2
echo gatk UnifiedGenotyper 2>>$log 1>&2
java -Xmx$java_memory -jar $gatk \
	-T UnifiedGenotyper \
	-R $ref_genome \
	-I 1.realn.recal.bam \
	-o 1.indel.vcf \
	-glm INDEL \
	-baq CALCULATE_AS_NECESSARY \
	--baqGapOpenPenalty 40.0 \
	--downsample_to_coverage 2000 \
	-nt 8 \
	-nct 2 \
	--min_indel_fraction_per_sample 0.25 \
	--min_indel_count_for_genotyping 20 \
	-stand_call_conf 50.0 \
	-stand_emit_conf 30.0 \
	-contamination 0.0005 \
	--min_base_quality_score 20 \
	--max_deletion_fraction 2 \
	-L ${intervals}


. $cmd_done
# BAQ gap open penalty # Phred-scaled gap open penalty for BAQ calculation. Although the default value is 40, a value of 30 may be better for whole genome call sets.

# max_deletion_fraction #If the fraction of reads with deletions spanning a locus is greater than this value, the site will not be considered callable and will be skipped. To disable the use of this parameter, set its value to >1.